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D0PEX-1D2C - A ONE-DIMENSIONAL, TWO -CONSTRAINT 
RADIATION SHIELD OPTIMIZATION CODE 
by Gerald P. Lahti 
Lewis Research Center 

SUMMARY 

A one -dimensional, two -constraint radiation shield weight optimization procedure 
and a computer program, DOPEX-1D2C, is described. DOPEX-1D2C uses the steepest - 
descent method to alter a set of initial (input) thicknesses of a spherical shield configura- 
tion to achieve a minimum weight while simultaneously satisfying two dose -rate con- 
straints. The code assumes an exponential dose-shield thickness relation with param- 
eters specified by the user. Code input instruction, a FORTRAN-IV listing, and a sam- 
ple problem are given. Typical computer time required to optimize a seven -layer shield 
is less than 1/2 minute on an IBM 7094. 


INTRODUCTION 

One -dimensional, single -constraint shield weight optimization procedures have been 
established for radiation shield design (refs. 1 and 2). Optimization procedures such as 
these take an initial configuration specified by the user and alter the thickness of each 
layer to achieve a minimum weight while simultaneously satisfying a single dose 
constraint. 

One method used to seek the minimum -weight configuration from a given starting 
configuration is that of steepest descent. Bernick (ref. 3) developed an algorithm to op- 
timize one -dimensional shields with one constraint by utilizing the steepest -descent 
method and programmed it as the OPEX code (ref. 4). OPEX was later reprogrammed 
for spherical geometry and reinterpreted as the OPEX-II code by Lahti (ref. 5). The 
steepest -descent strategy was extended in theory by Lahti (ref. 6) to the multidimen- 
sional, multiconstraint situation. A special case of this multidimensional optimization 
was programmed as the DOPEX code (ref. 7). DOPEX considered the special case 
where doses in each principal direction (i.e. , R-axis and plus and minus Z-axis of a 



right -circular cylinder) are dependent only on shield thicknesses in that direction. 
DOPEX might be considered to be a one -dimensional, one -constraint procedure coupled 
through the weight equation to operate in a two-dimensional, two- or three -constraint 
situation. 

The present report describes the next step in the evolution of multiconstraint opti- 
mization strategy, the DOPEX-1D2C code. DOPEX-1D2C is a one -dimensional, two- 
constraint optimization code and is presently coded for spherical geometry. The opti- 
mization algorithm and equations are those derived in reference 6, but because of an in- 
stability, a new algorithm for obtaining an initial configuration that would meet the input 
dose constraint had to be derived. Included in this report are derivations, data input in- 
structions, a FORTRAN -IV code listing, and a sample problem for DOPEX-1D2C. 


DOSE-THICKNESS RELATION 

‘ The DOPEX -1D2C code considers two dose rates D^(t) and DgCt) at some reference 
point in space. These dose rates are functions of layer thicknesses t. Each D m (t), 
m = 1 or 2, is in turn defined as 


IMAX 


m 


D m<« * £ D im<» 

i=l 


where 

^ i.U IL 

Di m (t) i component of m L total dose rate (e.g. , dose due to capture gammas from 
the first shield layer or dose due to inelastic gammas from the last shield 
layer, etc. ) 

D m (t) total dose rate, m = 1 or 2 (e.g. , might be total neutron dose rate and 
Dg might be total gamma dose rate) 

IMAX m number of dose -rate components in D m 

11 . 

t set of tj, the thickness of the j region 

Each dose -rate component D im (t) is further assumed to be of the form 


Dim* 4 ' * C im 



2 



where 


^im f^ed parameter 

N number of regions (thicknesses) in the configuration 

t. thickness of the j region, cm 

u. . ’’attenuation coefficient, " cm"^ 

^ijm 

This formulation is similar to that used in the OPEX-II and DOPEX codes and will 
not be detailed further in this report. For additional information on how to obtain these 
coefficients, the reader is referred to references 5 and 7. 


WEIGHT-MINIMIZATION PROCEDURE 

The procedure for obtaining the minimum -weight configuration by the method of 
steepest descent is presented in this section. Details of derivations can be found in ref- 
erence 6. 

The problem to be solved is that of minimizing the system weight w (a function of 
thickness t.) while simultaneously satisfying two dose -rate constraints Dj and D 2 (also 
a function of tj). Additional constraints to be imposed are that the thickness tj must be 
nonnegative and that certain thicknesses may be constrained to constant thickness. The 
problem, then, is to minimize w (tj, t 2 , . . . , t n ) with the following constraints: 

(1) Djd) = D° 

(2) D 2 (t) = D° 

(3) t j > 0 

(4) ^ is constant for any desired values of l . 

Constraint 3 ensures a physically meaningful solution. Constraint 4, fixed thick- 
ness, is useful if it is desired that some thicknesses be kept from changing during the 
course of the calculation (e. g. , the reactor core size and the reflector thickness). 

An n-dimensional Euclidean vector space with Cartesian coordinates, tj, t 2 , . . . , t^ 
is defined. The following vectors are defined on this space: 


3 




points in the direction of greatest weight decrease (steepest descent) along a hyperplane 
tangent to the hypersurface described by the equations 


Dj(t) = D° the first dose constraint 


apd 


Dg(t) = Dg the second dose constraint 


Components of u(t) represent increments of thickness to be added to t in order to ap- 
proach a minimum weight. 

The optimization algorithm, which is similar to that used in OPEX, OPEX n, and 
DOPEX, is as follows: 

(1) The unit vector u(t) is calculated for the present value of t. 

(2) A fraction, f, of u is added to t. The fraction f is an input parameter. (A 
value of f = 1. 0 has given satisfactory results. ) That is, 

4 



t = t + fu(t) 


(3) The new set of thicknesses t generally does not return the correct dose con- 
straints because of the nonlinear dose -thickness relation. A first-order correction dt 
is applied to t to restore the dose constraints. This correction (derived in appendix A 
as eq. (A 12)) is 

__ Jd^VDi • VP 2 ) - pj V 2 D 2 ] 7D 1 + K< VD 1 ' 7P 2 ) - 0*2 V 2 DllvD 2 

(vDj • vd 2 ) 2 - v 2 Dj v 2 d 2 

where is the difference, D m (t) - D^. Finally, then, in this step, t is incremented 
by dt to obtain a new t which satisfies dose constraints. 

(4) The system weight is calculated for this new t. If the change in weight (between 
that calculated in the present iteration and that in the past iteration) is less than some 
prescribed e, then the algorithm stops. Otherwise, steps 1 to 4 are repeated. 


INITIALIZATION PROCEDURE 

In general, an initial (or estimated) configuration will not satisfy dose constraints. 
Therefore, provision must be made to obtain some initial feasible solution before the op- 
timization can proceed. The method used in OPEX-II and DOPEX turned out to be un- 
stable for the present two -constraint case. Therefore, a new initialization algorithm 
was developed. Details of this derivation appear in Appendix A. 

The initialization algorithm proceeds as follows: 

(1) Calculate dose rates D^(t) and Dg(t) and compare with the required dose con- 
straints and Dg. 

(2) If D m (t) > 2D^ (where m = 1, 2), set = 0. 5 D m (t). 

If D m^ < °* 5 D m (where m = 1, 2), set = 2D m (t). 

U °-5 C [D m (t)/ D m] < 2 0 < where m = 1, 2), set = D m (t) - D° . 

(3) Using the differences Dj and Dg calculated in step 2, use equation (A12) to 

calculate dt and a new t = t + dt which will result in dose rates of approximately D. 

* 1 
and Dg . 
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(4) Calculate D m (t) and 


m 


(t) -D 


o 

m 



and compare to some convergence parameter e. If 


D (t) - D° 
m v m 


D 


m 


< e 


for m = 1 and 2, the algorithm is complete. If 




D ft) - D 
m' 


o 

m 



> e 


for m = 1 or 2, go to step 2. 

The present procedure has been uniformly convergent on all problems tested to date. 
The arbitrary range (factor of 0. 5 to 2.0), assumed for the first-order difference equa- 
tion to be valid, has so far been demonstrated to be adequate. The user may, of course, 
change this parameter in subroutine INIT if desired. 


D0PEX-1D2C CODE 

In this section, the details of the DOPEX-1D2C code are presented. Included are 
(1) an overview of the code, (2) a flow chart for data input, and (3) a sample problem and 

V 

sample problem output. A complete FORTRAN IV listing appears in appendix B. 




Overview of D0PEX-1D2C Code 

Objective . - DOPEX-1D2C is a radiation shield optimization code which obtains a 
minimum -weight configuration for a layered, spherical shell arrangement while simulta- 
neously satisfying two dose -rate constraints. This code alters layered configurations. 

It will not add shield layers, but it may delete layers. 

Method . - The method of steepest descent is used in an algorithm designed to alter 
an input configuration while simultaneously satisfying dose -rate constraints. 
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Required input data . - The user must provide an initial geometric configuration 
(layer thicknesses and densities of materials) and all data to be used in an empirical 
equation which relates changes in dose rates with changes in thickness. 

Computer and language used . - DOPEX-1D2C is written in ANSI-standard 
FORTRAN -IV and has been compiled and run on IBM 7094 and IBM 360/67 computers. 

Restrictions in problem size . - DOPEX-1D2C is presently set up for a maximum of 
25 regions and each of the two dose constraints is limited to 25 components. These 
limits could be increased as DOPEX-1D2C presently requires a small fraction of the 32K 
IBM 7094 available core. 

Typical machine time . - A typical DOPEX-1D2C run requires less than 1/2 minute 
on an IBM 7094. 

Precautions . - The final configuration predicted by DOPEX-1D2C is only as good as 
the empirical dose -thickness equation used in the calculation. The user is advised to 
check the predicted dose rate for the final configuration by some more exact method. 

Other geometries . - Other geometries can be used by rewriting subroutine WEIGHT 
to suit. 

More constraints . - DOPEX-1D2C is presently written to handle only two con- 
straints. However, equations for more constraints have been derived and could be in- 
corporated if deemed necessary. 


Flow Chart for Data Input 



The following is a flow chart for data input to DOPEX-1D2C. The symbol 
means read a card of set of cards as described. 




Card column 

Variable 

Description 

FORMAT 


1-80 


Title card (any alphanumeric 

20A4 

i 




information in card coIh 
umns .1 to 80) 


( 2 _ 



Control card 

215, 3E 10. 4 



1-5 

NREG 

Number of regions in prob- 






lem (<25) 




6-10 

MAX 

Maximum number of itera- 






tion steps allowed in opti- 
mization 




11-20 

EPS 

Convergence criterion for 



i 



weight (typically, 0.001) 
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Card column 

Variable 

Description 

FORMAT 

21-30 

EPSD 

Convergence criterion for 
initial dose (typically, 
0.001) 


31-40 

CON 

Fractional step size on u 
(typically, about 1.0) 



T(J) 

Thicknesses of each region; 
NREG values required 

7E10.4 


RHO(J) 

Densities of each region; 
NREG values required 

7E10.4 


NB(J) 

Constraint flags; NREG 
values required. 

1 ii. 

NB(J) = 0 constrains j 

region to constant thickness 

NB(J) = 1 allows j th region 
to change 

Dose rate control card 1 

2511 

15, E 10. 4 

1-5 

IMAX(l) 

Number of dose-rate com- 
ponents in first dose-rate 
constraint 


6-15 

DSTAR(l) 

Value of first dose-rate 
constraint 



EMU (I, J) 

Attenuation coefficients, 

For the first constraint, 
read a new card (or set of 
cards) with IMAX(l) values 
of p^, NREG cards or 
sets of cards are required. 

7E 10.4 


D(I, 1) 

Dose-rate components D. 
for first constraint; 
IMAX(l) values required 

7E10.4 
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Card column 


Variable 


Description 


FORMAT 



NIJ(I) Region number physically 1415 

associated with i in dose- 
rate component. This 
cross referencing is neces- 
sary to zero components if 
the j region is elinii- 
nated in optimization pro- 
cedure. A region may be 
identified more than once. 

IMAX(l) values required. 

Repeat cards 6 to 9 with 
data for second dose-rate 
constraint. 


Sample Problem 

A sample problem concerning a reactor, a molybdenum reflector, and a shield con- 
sisting of seven alternating layers of lithium hydride and tungsten is illustrated in fig- 
ure 1. This particular arrangement is the same as that used in the OPEX-II code and 
was optimized in that case for a single dose constraint of 2 mrem/hour at a distance of 
20 meters. 

As a test of the DOPEX-1D2C code, this same geometric layout is optimized for con- 
straints of 0. 2527 -mrem/hour neutron dose rate and 1. 74 73 -mrem/hour total gamma 
dose rate. These values were predicted by the single -constraint optimization calculation 
(OPEX-II). The present case, with two constraints, should result in a layout and weight 
similar to that predicted by the one -constraint optimization. The present case will be 
initiated at the same configuration as was the OPEX-II calculation. 


9 





The results of this sample problem are given in tables I to IV. Table I lists region 
descriptions, densities, initial thicknesses, and the thicknesses predicted by OPEX -II 
and by DOPEX-1D2C. Table n lists dose rates by component. Within the DOPEX-1D2C 
calculation, the neutron dose (item 1 in table II) is considered to be the first constraint 
and the total gamma dose (sum of items 2 to 12 in table n) is considered to be the second 
constraint. The configurations calculated by OPEX -II and DOPEX-1D2C are similar. 

7 

Although the thicknesses are not identical, the system weight, 3x10 grams, is. It has 
been observed that minimum -weight configurations for thick shields, as considered here- 
in, can exhibit considerable local variations. This simply demonstrates that the solution 
given by OPEX -II or DOPEX-1D2C is not a unique one. Table in lists the complete 
DOPEX-1D2C output for this case. Table IV gives the coefficients used in this run. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, May 3, 1973, 

503-05. 
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APPENDIX A 


RECALCULATION OF t TO MEET CONSTRAINTS 

The problem here is the following: given a vector t calculate 

D*(t) = Dj(t) - D° 

D*(t) = D 2 (t) - D° 

where D* (t) is the difference between the dose rate as calculated using t (i. e. , D(t)) and 

o * — 

the required dose constraint D . It is required that D (t) = 0 to satisfy dose -rate con- 

— * * 

straints. To obtain a feasible solution for t(i.e., one where Dj = 0 and Dg =0), as- 
sume the first-order correction is valid. That is, 

D* (! + dt) = D* (!) + VD* (!) • dt 


and 


D*(t + dt) =D*(t) + VD*(t) • dt 
For the present case we require dt for which the residuals 

D * (! + dt) = 0 

and 


D* (t + dt) = 0 

Also, to minimize perturbations on the system, we also will require that |dt| be a 
minimum. 

The problem is now restated as 

Minimize dt • dt (Al) 


subject to 


D*(t) + VD*(t) • dt = 0 


(A2) 
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and 


D*(t) + VD*(t) • dt = 0 

To solve this, form the Lagrangian <£ by taking 

if = dt • dt + r^D^t) + vd*( 7) • dt] + y 2 [ D 2^) + VD*(t) ■ dt] 


(A3) 


where y ^ and y 2 are multipliers to be determined. Equivalently, 

if = dt • dt + yj[Dj (t) + VD^tj • dt] + y^D* (7) + VD 2 (t) • dt] 


Take the derivatives of if with respect to the components of dt and set each equal to 
zero to obtain a stationary point: 


3jf 

aCdtp 


= 0 = 2 dtj 


+ y l VD 1 ' t l +y 2 VD 2 





5if 

a(dt 2 ) 


= 0 = 2 dtg + y j VDj 


' ‘ 2 - V 2 VD 2 





(A4) 


i |^ = 0=2<« N *r 1 VD 1 -{ N+y2 VD 2 -‘t N 


J 


A f Jl — 

where tj are unit vectors in the i direction. The arguments, t, on Dj and D 2 will 
presently be dropped. 

Equation (A4) can be written as 


3if 

3(dt) 


= 0 = 2 dt + y^ VDj + y 2 VDg 


(A5) 


Equations (A2) to (A4) form a set of N + 2 equations in N + 2 unknowns (dt. and y ^ 
and y 2 ). To solve this, from equation (A5) obtain 


1 

2 


dt = - — (y x VDj + y 2 VD 2 ) 


(A6) 
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Substitute equation (A6) into (A2) and (A3) and obtain 


* vd i ■ (• zK vd i + y 2 ^ = ° 

and 

D 2 t7D 2-(4) <y l VD l +,, 2 VD 2 )=0 

From equation (A8), 

. 2^2 - >-2 V \ 

Yl ~ VD V ^2 

Substitute equation (A9) into (A7) and solve for to obtain 

2D*(VDj • VD 2 ) - 2D* V 2 Dj 

2 ( VD 1 ' ™2f - v2 °l V \ 

Substitute equation (A 10) into (A9) to obtain 

2D*(VDj • VD 2 ) - 2D* V 2 D 2 
1 ( VD 1 ' VD 2) 2 - v2 °l v2 °2 

Finally, use equations (A6), (A 10), and (All) to obtain the value of dt: 

[ D 2<VDr ^P 2 l- P t v2p 2 1 VP l - [ D 1< VP 1 * vp 2)- p 2 ?2 d i] VD 2 

(VD, • VD 2 ) 2 - V 2 D, V 2 D 2 


(A7) 


(A 8) 


(A9) 


(A 10) 


(All) 


(A 12) 
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APPENDIX B 


FORTRAN IV LISTING 


The FORTRAN -IV listing for the DOPEX-1D2C code is presented in this appendix. 
Included are the main routine and subroutines INPUT, INIT, CLEAR, DOZE, and 
WEIGHT. The coding is ANSI standard FORTRAN -IV and has been run successfully on 
both the IBM 7094-11 and IBM 360/67 computers. 


Description of Subroutines 

The subroutines and their function in DOPEX-1D2C are as follows: 

(1) Main program: The main program calculates the steepest -descent vector and 
alters shield thicknesses; it calls subroutines INPUT, INIT, DOZE, WEIGHT, and 
CLEAR. 

(2) Subroutine INPUT: This subroutine reads all input data. 

(3) Subroutine INIT: If initial thicknesses do not result in the dose rates required by 
the constraints, INIT is called to calculate a set of thicknesses which does, so that the 
optimization can begin. This calculation is carried out without regard to weight minimi- 
zation. INIT calls subroutines DOZE, WEIGHT, and CLEAR. 

(4) Subroutine WEIGHT: This subroutine calculates the total weight and the deriva- 
tives of weight with respect to thickness. It is presently coded for spherical geometry. 

(5) Subroutine DOZE: This subroutine calculates dose rates, the derivatives of dose 
rates with respect to thickness, and products of dose and weight derivatives. 

(6) Subroutine CLEAR: If a shield layer is eliminated in the course of a calculation, 
subroutine CLEAR is called to set that layer thickness to zero for the remainder of the 
calculation, to set the constraint flag so that it is maintained at zero thickness, and to 
zero out a dose component if it originated in the removed layer. 


Program Listing 


C NEW ONE-DIMENSIONAL, DOUBLE CONSTRAINT OPTIMIZATION CODE 

C ...GPL 3 JAN 1972... REVISED DEC 1972... 

COMMON MAX, EPS, EPSD, CON, C A, WT, GO, NREG, 

1 DSTAR( 2 ) , T ( 25 ) , RHO(25), NB(25), IMAX(2), 0(25,2), 

2 EMU (2 5, 2 5, 2 ) , NIJ(25,2), D0SE(2), C(25,2), A(25,2), 

3 AA ( 2 ) , AXA, AG (2 ) , U(25), 0(25) 

CCC 

C MAX NUMBER OF ITERATIONS ALLOWED 

C EPS PRECISION REQD ON WE I GHT CONVERGENCE 
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o o o o a o a o o o o o o o o o o o o o o o d n o r> o 


EPSD 

CON 

CA 

WT 

DSTAR(2) 
NREG 
T(25 ) 
RHO(25 ) 

NB (25 ) 

I MAX (2) 
D(25/2) 
EMU(25,25, 
N IJ (25,2) 
HOSE (2 ) 

C (25* 2 ) 
A(25/2) 

AA (2 ) 

AXA 
AG (2) 

U( 25 ) 

G ( 25 ) 


PRECISION ON DOSE IN INITIALIZATION 

FRACTION FOR OPTIMIZATION UNIT VECTOR 

NOT USEO PRESENTLY 

TOTAL SYSTEM WEIGHT 

G .DOT. G 

DOSE CONSTRAINTS 

NUMBER OF REGIONS 

THICKNESS OF EACH REGION 

DENSITY OF EACH REGION 

CONSTRAINT FLAG FOR EACH REGION 

NUMBER OF DOSE COMPONENTS IN EACH DOSE CONSTRAINT 
VALUE OF THE RESPECTIVE DOSE COMPONENTS 
2) THE COEFFICIENTS FOR EACH REGION AND EACH CONSTRAINT 
CROSS-REFERENCE-REGION NUMBER FOR EACH DOSE COMPONENT 
PRESENT EVALUATION OF THE DOSE FOR EACH CONSTRAINT 
LEADING COEFFICIENT IN 0(25, 2) EQUATION 
D-D( I )/D-T(J) . . . 1 = 1*2/ J = !t/ NREG 
DEL-SQUARED OF D ( I ) 

DEL-D(l) .DOT. DEL-D(2) 

DEL-D(I) .DOT. DEL-W 

COMPONENTS OF OPTIMIZATION UNIT VECTOR 
D-W/D-T (J ) 


5 CALL INPUT 
CALL WEIGHT 
CALL DOZE 
WR I TE (6, 7 ) WT 

7 FORMAT (//20H0*** SYSTEM WEIGHT =/lPE12.4) 

IF((DSTAR(1) .GT. 0.0) .AND. (DSTAR ( 2 ) .GT. 0.0)) GO TO 10 
DSTAR(l)=DOSE(l) 

DSTAR ( 2 ) =DOSE ( 2 ) 

WRITE (6/ 8) (DSTAR( 1C), IC=1,2) 

8 FORMAT (44H0*** PROBLEM WILL OPTIMIZE TO CONSTRAINTS OF ,1P2E12.4) 
GO TO 12 


10 CALL INIT 
WR I TE (6/ 7) WT 

12 WR I TE (6, 15 ) 

15 FORMAT (40H1******** BEGIN OPTIMIZATION *********** ) 

MAIN PROGRAMME 

DO 100 K=1/MAX 
WTX=WT 

DET=AA ( 1 ) *AA ( 2 ) -AXA* *2 
ALF1= (AG (2 ) * AX A- AG (1)*AA(2))/DET 
ALF2= ( AG (1 )* AX A- AG (2)*AA(1))/DET 

F0RB2=GG+ALF1* ( 2 . 0*AG (1)+ALF1*AA(1)+2.0*ALF2*AXA) 
FORB2=FORB2+ALF2*(2.0*AG(2)+ALF2*AA(2)) 

C THE ABOVE IS EQUATION(20) FROM NASA TM X-2435 4*BETA**2 

DENOM=SQRT(FORR2) 

C 

DO 20 I =1/NREG 

I F (NB ( I ) .EQ. 0) GO TO 20 

U ( I ) = (G ( I)+ALF1*A(I/1)+ALF2*A( 1,2) )/DENOM 

T(l)=T(l)-CON*U(l) 

I F (T ( I ) .LE. 0.0) CALL CLEAR(I) 

20 CONTINUE 
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RESTORE ROSES TO CONSTRAINT 
CALL OOZE 

DENOM=AXA**2-AA(l)*AA(2) 

D1=D0SE(1)-DSTAR(1) 

D2=D0SE ( 2)-DSTAR(2) 

GAMMA1=(D2*AXA-D1*AA(2) J/DENOM 

GAMMA2=(D1*AXA-D2*AA(1) )/DEN0M 

DO 30 I =1/NREG 

IFCNB(I) .EQ. 0) GO TO 30 

T ( I ) =T ( I ) -GAMMA 1*A (1,1) -GAMMA 2*A( I ,2) 

I F (T ( I ) .LE. 0.0) CALL CLEAR(I) 

30 CONTINUE 

CALL WEIGHT 
CALL DOZE 

WRITE(6,40) K, WT, (DOSE(IC), 10=1,2) 

40 FORMAT (8H0*** K =,13, 5X, 12H*** WEIGHT =,1PE12.4, 

1 5X, 17H*** TOTAL DOSES =,2E12.4) 

WRITE(6,42)(T(I), l=l,NREG) 

42 FORMATOH T = , 1P10E12 . 4/ (9X, 10E12 . 4 ) ) 

I F ( ( K/ 5 ) *5 .NE. K) GO TO 55 
ISKI P=1 

45 DO 48 IC-1,2 
I MMX= I MAX (1C) 

48 WR I TE (6, 49 ) 1C, (0(1,10,1 = 1, I MMX ) 

49 FORMAT (42H0 DOSE COMPONENTS FOR CONSTRAINT NUMBER, 12/ 

1 (9X, 1P10E12 . 4 ) ) 

GO TO(55,5), ISKIP 

55 I F ( ABS( (WTX-WT ) /liTD .GT. EPS) GO TO 100 
I SK I P=2 
GO TO 45 

100 CONTINUE 
STOP 
END 


SUBROUTINE INPUT 

NEW ONE-DIMENSIONAL, DOUBLE CONSTRAINT OPTIMIZATION CODE 
...GPL 3 JAN 1972... REVISED DEC 1972.. 

COMMON MAX, EPS, EPSD, CON, CA, WT, GG, NREG, 

1 DSTARC2), T(25 ), RHO(25), NB(25), I MAX ( 2 ) , D(25,2), 

2 EMU (25, 25 , 2 ) , NIJ(25,2), D0SE(2), C(25,2), A(25,2), 

3 AA ( 2 ) , AXA, AG (2 ) , U(25), G(25) 


MAX 

EPS 

EPSD 

CON 

1C 

DSTAR(2) 
NREG 
T ( 25 ) 
RH0( 25 ) 


NUMBER OF ITERATIONS ALLOWED 
PRECISION REQD ON WEIGHT CONVERGENCE 
PRECISION ON DOSE IN INITIALIZATION 
FRACTION FOR OPTIMIZATION UNIT VECTOR 
CONSTRAINT COUNTER 
DOSE CONSTRAINTS 
NUMBER OF REGIONS 
THICKNESS OF EACH REGION 
DENSITY OF EACH REGION 
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NB(25) CONSTRAINT FLAG FOR EACH REGION 

I MAX (2) NUMBER OF DOSE COMPONENTS IN EACH DOSE CONSTRAINT 
0(25,2) VALUE OF THE RESPECTIVE DOSE COMPONENTS 
EMU(25, 25,2) THE COEFFICIENTS FOR EACH REGION AND EACH CONSTRAINT 
NIJ(25,2) CROSS-REFERENCE-REGION NUMBER FOR EACH DOSE COMPONENT 
D0SE(2) PRESENT EVALUATION OF THE DOSE FOR EACH CONSTRAINT 

C ( 2 5 , 2 ) LEADING COEFFICIENT IN D(25,2) EQUATION 

DIMENSION ALPHA ( 20 ) 

READ (5, 5 ) ALPHA 
5 FORMAT(20A4) 

WRlTE(6,10) ALPHA 
10 FORMAT (1H1,5X,20A4) 

READ(5, 15) NREG, MAX, EPS, EPSD, CON 
15 FORMAT (215, 4E10.4) 

WRITE(6,20) NREG, MAX, EPS, EPSD, CON 
20 FORMAT ( 8H0 NREG = , I 3, 5X, 5HMAX = , 1 4 , 5 X , 5HEPS =,F8.5,5X, 

1 6HEPSD =,F8.5,5X, 5HCON =,F8.3,5X) 

C 

READ(5, 30) (T ( I ) , 1-1, NREG) 

30 FORMAT (7E10 . 0 ) 

READ(5, 30) (RHO(I), 1=1, NREG) 

READ(5, 35 ) (NB(I), 1=1, NREG) 

35 FORMAT (2511) 

WR I TE (6, 36 ) ( I , T( I ) , RHO(I), NB(I), 1=1, NREG) 

36 FORMAT (//34H0REG I ON T( I ) RHO(I) NB ( I ) / ( I 7 , 2F10 . 3, I 7 ) ) 

C 

DO 70 I C = l, 2 

READ(5, 37) IMAX(IC), DSTAR(IC) 

37 FORMAT ( I 5, E10 . 4 ) 

WR I TE (6, 38 ) 1C, IMAX(IC), DSTAR(IC) 

38 FORMAT (9H1****IC =, I 2, 5X, 6H I MAX =, I 3, 5X, 12HCONSTRA I NT =,1PE12.4) 

I MMX= I MAX (1C) 

C 

WRI TE (6, 39 ) 

39 FORMAT (27H0REG I ON- J MU(I,J)) 

DO 40 J=l, NREG 

' READ (5, 30 ) (EMU( l,J, 1C), 1=1, IMMX ) 

40 WRITE(6,45) J, (EMU( I , J, 1C) , 1=1, IMMX) 

45 FORMAT (16, 3X, 1P9E12 . 4/ (9X, 1P9E12 . 4 ) ) 

C 

READ (5, 30 )(D(I,IC), 1=1, IMMX) 

RE AD (5, 50 )(NIJ( I, 1C), 1=1, IMMX) 

50 FORMAT (14 I 5 ) 

C 

DOSE ( 1C ) = 0 . 0 

DO 60 1=1, IMMX 

DOSE ( I C )=DOSE ( I C )+D( I , 1C) 

BB=0 . 0 

DO 55 J=l, NREG 
55 BB=BB+EMU( I , J, IC)*T(J) 

60 C(l, IC)=D(I, IC)*EXP(BB) 

WRITE (6,65) (l,C( I, IC),D(I, I C ) , N I J ( 1 , 1C), 1 = 1, IMMX) 
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65 FORMAT (1H0/34H0 I C(l) 0(1) N I J/ ( I 5, 1P2E12 . 4, I 5 ) ) 

WRITE (6,68) DOSE (1C) 

68 FORMAT ( 16H0 TOTAL OOSE =,1PE12.4) 

C 

70 CONTINUE 
100 RETURN 
END 


SUBROUTINE IN IT 

C \ NEW ONE-DIMENSIONAL, DOUBLE CONSTRAINT OPTIMIZATION CODE 
C . . .GPL 3 JAN 1972. . . 

C ...REVISED GPL 18 DEC 1972... 

COMMON MAX, EPS, EPSD, CON, CA, WT, GO, NREG, 

1 DSTAR( 2 ) , T ( 25 ) , RHO(25), NB(25), IMAX(2>, D(25,2), 

2 EMU(25, 25,2), NIJ(25,2), OOSE(2), C(25,2), A(25,2), 

3 AA(2) , AXA, AG ( 2 ) , U(25), G(25) 

ccc 

DIMENSION DEL ( 2 ) 

C 

1 FACTOR=0 . 5 
FLOW=FACTOR 
FHI=1.0/FACTOR 
C 
C 

DO 67 K=l, MAX 
C 

DO 55 I C = l, 2 

IF(DOSE(IC)-FHI*DSTAR(IC)) 52,52,51 

51 DEL ( I C )=DOSE ( IC)*(1.0-FLOW) 

C DOSE (1C) IS TOO HIGH. . .ADJUST T CORRESPONDING TO 
C A DOSE CONSTRAINT OF FLOW*OOSE( «C) 

GO TO 55 
C 

52 I F (DOSE ( I C ) -FLOW*DSTAR( I C ) ) 53,53,54 

53 DEL ( I C ) =DOSE (IC)*(1.0-FHl ) 

C DOSE IS TOO LOW... ADJUST T CORRESPONDING TO 

C A DOSE CONSTRAINT OF FHI*DOSF.(IC) 

GO TO 55 

54 DEL ( I C ) =DOSE ( I C )-DSTAR( I C ) 

C DOSE( IC)/DSTAR( |C) IS .GT. FACTOR .AND. .LT. 1/FACTOR 

C ASSUME FIRST ORDER CORRECTION MODEL TO WORK 

55 CONTINUE 
C 

DENOM=AX A**2-AA ( 1 ) *AA ( 2 ) 

GAMMA1= (DEL ( 2 ) *AXA-DF.L (1 ) *AA (2 ) )/DENOM 
GAMMA2= (DEL ( 1 ) *AXA-OEL (2 ) *AA(l))/DENOM 
C 

U ( I )=GAMMA1*A (1,1 )+GAMMA2*A ( 1,2) 

C 

DO 66 1=1, NREG 

I F ( M B ( I ) .EQ. 0) GO TO 66 

U( I )=GAMMA1*A( I , 1 )+GAMMA2*A( 1,2) 

T( I )=T ( I )-U( I ) 

I F (T ( I ) .LE. 0.0) CALL CLEAR(I) 

66 CONTINUE 
C 

C TEST CONVERGENCE 
CALL DOZE 
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I F (ABS (DOSE (1 )/DSTAR Cl ) -1.0) .OT. EPSD) GO TO 6? 
I F (ABS (DOSE ( 2 )/DSTAR (2 ) -1.0) .LT. EPSD) 00 TO 60 
67 CONTINUE 
K=777 


UPDATE ALL DERIVATIVES 

69 CALL WEIOHT 
CALL DOZE 

WRI TE (6, 70 ) K, WT, (DOSE(IC), IC-1,2) 

WRITE (6/75) (TO), l=l,NREG) 

DO 80 IC=l/2 
I MMX 3 I MAX (1C) 

80 WRITE (6, 85) 10/(0(1,10/ 1-1, IMMX) 

RETURN 

70 FORMAT (28H1*** IN ITIALI ZAT I ON REQUIRED , I3,12H ITERATIONS./ 

1 12H0 WEIOHT =, 1PE12 .4/ 

2 12H0 DOSES ARE /1P2E12.4) 

75 FORMAT ( 3 2H0**** INITIAL THICKNESSES ARE . . . / ( 1P10E12 . 4 ) ) 

85 FORMAT (41H0**** INITIAL DOSES FOR CONSTRAINT NUMBER, 12/ 

1 ( 1P10E12 . 4 ) ) 

END 


SUBROUTINE OOZE 

NEW ONE-DIMENSIONAL, DOUBLE CONSTRAINT OPTIMIZATION CODE 
...GPL 3 JAN 1972... 

COMMON MAX, EPS, EPSD, CON, CA, WT, GO, NREG, 

1 DSTAR(2), T (25 ) , RHO(25), NB(25), IMAX(2), D(25,2), 

2 EMU(25, 25,2 ), NIJ(25,2), D0SE(2), C(25,2), A(25,2), 

3 AA(2) , AXA, AO ( 2 ) , U(25), 0(25) 

ccc 

DO 30 IC-1,2 
DOSE( 10 = 0.0 
I MMX= I MAX (1C) 

DO 10 1=1, IMMX 
BB=0 . 0 

DO 5 J=1,NRE0 
5 BB=BB+EMU( I , J, IC)*T(J) 

D(l, I C ) =C ( I , IC)*EXP(-BB) 

10 DOSE ( I O-DOSE ( I C)+D( 1 , 1C) 

C 

AA( 10 = 0.0 
AG( 10 = 0.0 
DO 20 K=1,NRE0 
A(K, IC)=0.0 

I F (NB (K) .EQ. 0) 00 TO 20 

DO 15 1=1, IMMX 

15 A (K, I C )-A (K, IC)-EMU( I , K, IC)*D( I, 1C) 

AA ( I C )=AA ( I C)+A (K, 10**2 

C *** 

C *** MAKE SURE THAT YOU HAVE CURRENT VALUES OF AO AND 0 

C *** WHEN EVALUATING THIS NEXT TERM 

AG( IC)=AG(IC)+A(K, I C) *0 (K) 

20 CONTINUE 
C 
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30 CONTINUE 
C 

AXA=0 . 0 

DO 40 K=1,NP.EG 
40 AXA=AXA+A(K,1)*A(K,2) 
C 

RETURN 

ENO 


SUBROUTINE WEIGHT 

NEW ONE-DIMENSIONAL, DOUBLE CONSTRAINT OPTIMIZATION CODE 
...GPL 3 JAN 1972... 

COMMON MAX, EPS, EPSD, CON, CA, WT, GO, NREG, 

1 DSTAR( 2 ) , T ( 25 ) , RHO(25), NB(25), IMAX(2), D(25,2), 

2 EMU( 25/25,2), NIJ(25,2), DOSE(2), C(25,2), A(25,2), 

3 AAC2), AXA, AG(2), U(25), G(25) 

DATA FORPI/12. 56636/ 

SPHERICAL GEOMETRY 
DIMENSION R(26) 

R( 1 ) =T ( 1 ) 

RRR=R( 1 ) **3 
WT=RHO(l)*RRR 
DO 5 1=2, NREG 
R3=RRR 

R( I ) = R ( I -1 ) +T ( I ) 

RRR=R ( I ) **3 
WT=WT+RHO( I ) * ( RRR-R3 ) 

5 CONTINUE 
WT=WT*FORP 1/3.0 

CALCULATE PARTIAL WEIGHT DERIVATIVES 
GG = 0 . 0 

DO 7 1=1, NREG 
RR=0 . 0 
G ( I )=0.0 

I F (NB ( I ) .EQ. 0) GO TO 7 
DO 6 J= I ,NR 
R2 = RR 

RR=R(J )**2 

6 G( I )=G( I )+RHO(J)*(9R-R2) 

G ( I ) =G ( I ) *FORP 1 
GG=GG+G( I )**2 

7 CONTINUE 
RETURN 
END 
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SUBROUTINE CLEAR(J) 

NEW ONE-DIMENSIONAL, DOUBLE CONSTRAINT OPTIMIZATION CODE 
...GPL 3 JAN 1972 .. . 

COMMON MAX, EPS, EPSD, CON, C A, WT, 00, NREO, 

1 DSTAR(2 ) , T( 25 ) , RHO(25), NB(25), I MAX C 2 ) , D(25,2), 

2 EMU (25, 25 , 2 ) , N!J(25,2), D0SE(2), C(25,2), A(25,2), 

3 AA (2 ) , AXA, AG( 2 ) , U(25), 0(25) 


THE JTH REGION HAS JUST BEEN WIPED OUT 


T(J )=0 . 0 
NB (J )=0 
DO 5 I C=l, 2 
I M- I MAX ( I C ) 

DO 5 1=1/ I M 

I F ( J .EQ. N I J ( I , I C ) ) 0(1,10=0.0 

5 CONTINUE 
RETURN 
END 
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TABLE I. - REGION DESCRIPTION 


Region, 

j 

Description 

Density , 
g/cm 3 

Initial thickness, 
cm 

(guess) 

Final prediction of 
shield thicknesses 
by OPEX-H, 
cm 

Final prediction of 
shield thicknesses 
by DOPEX-1D2C, 
cm 

1 

Reactor core 

9.957 

26. 0 radius 



2 

Plenum 

8.647 

2. 50 



3 

Pressure vessel 

16.763 

.60 



4 

Molybdenum reflector 

9.234 

11.00 



5 

Lithium hydride 

.75 

(17.90) 

20.52 

22.49 

6 

Tungsten 

19.3 

(7. 00) 

9.71 

9.62 

7 

Lithium hydride 

.75 

(14.00) 

12.32 

12.75 

8 

Tungsten 

19.3 

(5.00) 

2.82 

2.99 

9 

Lithium hydride 

.75 

(10.00) 

10.32 

9.51 

10 

Tungsten 

19.3 

(3. 50) 

2.33 

1.81 

11 

Lithium hydride 

.75 

(59.50) 

39.29 

38.24 


TABLE H. - DOSE COMPONENTS 


Item, 

i 

Dose component, 
D i 

Initial shield 
thicknesses 

Final shield thicknesses 

OPEX-n 

DOPEX-1D2C 

Value of dose component D., 
mrem/hr 

i 

Neutron 

0.02430 

0.2527 

0.2527 

2 

Core gamma 

. 00303 

.0079 

.0110 

3 

Plenum, pressure vessel capture gamma 

.00196 

.0049 

.0069 

4 

Plenum, pressure vessel inelastic gamma 

. 00220 

.0055 

.0077 

5 

Reflector capture gamma 

.204 

.4782 

.6710 

6 

Reflector inelastic gamma 

. 00504 

.0129 

.0178 

7 

Region 6 tungsten capture gamma 

.0921 

.4007 

.2703 

8 

Region 6 tungsten inelastic gamma 

. 00974 

.0949 

.0842 

9 

Region 8 tungsten capture gamma 

.0988 

.2696 

.2590 

10 

Region 8 tungsten inelastic gamma 

.0278 

.0759 

.0788 

11 

Region 10 tungsten capture gamma 

.201 

.2650 

.2332 

12 

Region 10 tungsten inelastic gamma 

.0947 

. 1320 

. 1073 

Total 

0.7647 

2. 000 

2.000 
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TABLE HI. - SAMPLE D0PEX-ID2C OUTPUT 
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TOTAL DOSE = 2.4300E-02 



IMAX = 11 CONSTRAINT = 1.7473E 00 
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DOSE COMPONENTS FOR CONSTRAINT NUMBER 2 

I.1021E-02 5.8R93E-03 7.7329E-03 6.7103E-01 1.78165-02 2.70355-01 8.<323AE-02 2.5932E-31 7.979SE-02 2.3318E-01 

l. 0725E-01 
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Figure 1. - Geometry for sample problem. 
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